
dxf = 1.0*sqrt(udata.dx(1)*udata.dx(1)+udata.dx(2)*udata.dx(2));
udata.N_fractures = size(frac_info,1);

lamda=0.95;        % decay rate of fracture size if oversize by jiang 2022-5-11
kesi=0.01;         % a low angle to avoid vertical and horizontal fracture with zero dxi and dyi by jiang 2022-5-11

xc = []; yc = [];
xb = []; yb = [];
xe = []; ye = [];

xy=[frac_info(:,1),frac_info(:,3)];                                         % start points of each lien
beta=90-frac_info(:,6);                                                     % Fracture strike 90 in north; 0 for east, -90 for south 
L=round(frac_info(:,5));                                                    % length of each fracture line
NL = floor(L./dxf);                                                         % number of segement 

for i=1:udata.N_fractures 
    pass = 0;
    ra=1.0;
    r=NL(i);
    while pass < 1
        
        x = xy(i,1);
        y = xy(i,2);
        r = floor(r*ra);
        
        dyi = sin((beta(i)+kesi)*pi/180)*dxf;                                % jiang 2022-5-11  size in y direciton of each segment
        dxi = sqrt(dxf*dxf-dyi*dyi);                                         % size ub x direction of each segment

        xbi = x:dxi:x+r*dxi-dxi;                                             % coordiantes of start x of each segment
        xei = x+dxi:dxi:x+r*dxi;                                             % coordinates of end x of each segment
        ybi = y:dyi:y+r*dyi-dyi;                                             % coordinates of start y of each segment
        yei = y+dyi:dyi:y+r*dyi;                                             % coordinates of end y of each segment

        if (all(xbi <= udata.len(1)) && all(xei <= udata.len(1)) && ...
            all(ybi <= udata.len(2)) && all(yei <= udata.len(2)) && ...
            all(xbi >= 0) && all(xei >= 0) && ...
            all(ybi >= 0) && all(yei >= 0))
            pass = 1;                                                        % all coordinates in the model domain
        end
        ra=lamda;                                                            % jiang 2022-5-11
    end
    
    udata.Nf_i(i) = length(xei);                                             % segments of each fractures
    
    xb =[xb round(xbi .* 1e4) ./ 1e4];                                       % keep four digits after point
    xe =[xe round(xei .* 1e4) ./ 1e4];
    yb =[yb round(ybi .* 1e4) ./ 1e4];
    ye =[ye round(yei .* 1e4) ./ 1e4];
end

udata.Nf_f = length(xe);                                                     % total number of fracture segment

XY1 = [xb' yb' xe' ye'];
XE = [xe; ye]';

udata.frac_angle = atan2((ye-yb),(xe-xb))*180/pi;
udata.dxf = dxf;